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Abstract 

The convergence problem of the Laplace-Beltrami operators plays an 
essential role in the convergence analysis of the numerical simulations of 
some important geometric partial differential equations which involve the 
operator. In this note we present a new effective and convergent algorithm 
to compute discrete Laplace-Beltrami operators acting on functions over 
surfaces. We prove a convergence theorem for our discretization. To 
our knowledge, this is the first convergent algorithm of discrete Laplace- 
Beltrami operators over surfaces for functions on general surfaces. Our 
algorithm is conceptually simple and easy to compute. Indeed, the con- 
vergence rate of our new algorithm of discrete Laplace-Beltrami operators 
over surfaces is 0(r) where r represents the size of the mesh of discretiza- 
tion of the surface. 

Keywords:Local tangential polygon; discrete Laplace-Beltrami operators; 
Configuration Equation 



1 Introduction 

Let S be a smooth surface in the 3D space. The Laplace-Betrami (LB) operator 
is a natural generalization of the classical Laplacian A from the Euclidean space 
to S. It is well-known that the LB operator is closed related to the mean 
curvature normal by the relation As(p) = 2H{p). The LB operator plays 
important role not just in the study of geometric properties of E, but also in 
the investigation of physical problems, like heat flow and wave equations, on 
E. Moreover, the LB operator has recently many applications in a variety of 
different areas, such as surface processing [6l[T3], signal processing [Ml [T5l [TBI 
[T5] and geometric partial differential equations [T1[TT1[T2]. Since the objective 
underlying surfaces to be considered are usually represented as discrete meshes 
in these applications, there are tremendous needs in practice to discretize the 
LB operators. 

Even though the computation of the LB operators is important for many 
applications, there does not exist a simple "convergent" discrete approximation 
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of the LB operators for general surfaces. In this paper we shall present a new 
effective and convergent algorithm to compute discrete Laplace-Beltrami oper- 
ators acting on functions over surfaces. In fact, we shall prove the following 
convergence theorem. 

Main Theorem. Given a smooth function h on a regular surface S and a 
triangular surface mesh S = (V, F) with mesh size r, one has 

/^^h{v) ^ /^Ah{v) + 0{r) (1) 

where the discrete LB operator A^/i(u) is given in Equation i2^) . 

We shall give a mathematical proof of this convergence result. To our knowl- 
edge, this is the first convergent algorithm of discrete Laplace-Beltrami opera- 
tors over surfaces for functions on general surfaces. The idea of our algorithm 
can be divided into two parts: First, we shall introduce a notation of the local 
tangential polygon and lift functions and vectors on a triangular mesh, obtained 
from the discretization of the surface under consideration, to the local tangential 
polygon, and second, we shall give a new method to define the discrete Laplace- 
Beltrami(LB) operator acting on functions on a 2D polygon. Our algorithm 
is conceptually simple and easy to compute. The convergence rate of our new 
algorithm of discrete Laplace-Beltrami operators over surfaces is 0{r) where r 
represents the size of the mesh of discretization of the surface. We also present 
our numerical results to support this in section 4. 



2 Laplace-Beltrami operator and its discretiza- 
tions 

Let S be a regular surface in the 3D Euclidean space R^. Consider a local 
parameterization /i : J7 — S with h{ui,U2) = (a;(iti, M2), M2), U2)) € 
S, (wi,M2) e C/ C M^. For the details, we refer to do Carmo [HI IS]- Then the 
Laplace-Beltrami operator As applying to a function / on E is given by 



d 

dui 



df 

dj 



where (.9''') is the inverse of the matrix {gij) with g 

dh dh 



det {gij) and 



gij 



dui ' du. 



> 



(2) 



(3) 



Consider a triangular discretization S — {V,F) of the surface S, where 
V — {vi\l < i < ny} is the list of vertices and F = {Tfe|l < k < np} is the list 
of triangles. Let w be a vertex in V and N{v) the index of one-ring neighbors 
of the vertex v . Next we recall several discretizations of A/ for a function 
/ on as follows. For more discussions, see also Xu |19l [20] . 



2.1 Taubin's et al. Discretization 

Taubin considered in |16J the following form of discretization of A/: 

i£N(v) 



(4) 
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ieN{v) 



— 1. There are 



where the weights Wi are nonnegative numbers with ^ 

several choices for the weights Wi. An obvious choice is the uniform weights 
where |-/V(w)| is the cardinahty of the set N{v). A general way to 



determine the weights uji is to use the following formulation: 



UJi 



'Ek€N(v) (t>{v,Vk) 



with a nonnegative function (j){v,Vi). Fujiwara takes (j)(v,Vi) = 
brun's et al. 7 defines the weights cji as 



cot a.i 



- cot /3i 



'keN(v 



cot a/c + cot Pk 



where and /3i are the triangles as shown Figure [TJ 



(5) 



Des- 



(6) 




Figure 1: The angles and /3i. 

It is obvious that the discretization of A/ can not be a correct approx- 
imation of A/ since it approaches zero as the size of the surface mesh goes to 
zero. 

2.2 Mayer's et al. Discretization 

For a function / on S, Green's formula gives 



Af{x)dx= / dnf{s)ds (7) 

D(z,e) JdD(z,t) 

where D{z^ e) is a small disk at a point z on the surface E, and n is the 
intrinsic outer normal of the boundary of the disk. Mayer discretized ([7]) at w 
over the triangular surface mesh S and obtained the following approximation 

where A{v) is the sum of areas of triangles around v, and k,m Cz N{v) D N{vi). 
It can be checked directly that the formula (O is derived from ([7]) by ap- 
proximating J^i^^^^^Asf{x)dx, dnf{s) and ds with Af{v)A{v), and 

— — }hlJA!hs.^z]hl ^ respectively. Therefore, the discretization in ([5]) is an approx- 
imation of As/ at V. 
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2.3 Desbrun's et al. Discretization 



It is well-known in the theory of differential geometry that the mean curvature 
normal satisfies the following formula: 

lim ^ = -H{p) (9) 

where A is the area of a small region around the point p, and V is the gradient 
with respect to the (x, y, z) coordinates of p. From Equation ([S]), Desbrun et al. 
got the following approximation: 

^/w - E - /(.)|. (10) 

^ ' ieN{v) 

where N{v) is the index set of 1-ring neighboring vertices of vertex v, ai and Pi 
are as in ^ and A{v) is the sum of areas of triangles around v. 



2.4 Xu's Discretization 

In 2004, Xu presented two discrete Laplace-Beltrami methods in a triangular 
mesh from Green's formula and the quadratic fitting. Following Equation ([7]), 
xu introduced his discretization 

A/(^;) = ^ E < V/(p) + V/(K),^^>lk^-«.JI (11) 

where V/(t;) is the gradient of / at v, A{v) is the sum of area of the triangles 
that contain v and i>i is the unit outward normal of the edge ViVi_^ . Xu also 
used the biquadratic fitting of the surface data and function data to calculate 
the approximate LB operator. He introduced complexity weights of the equation 
(jH). This kind of weights can be found in [20J. 

The convergence problem of these discrete LB operators over triangular sur- 
face meshes has been investigated by Liu, Xu and Zhang in [Tni HH [10] ■ None of 
the above mentioned discretizations of the LB operators has ever been proved 
to be convergent for general surfaces and functions. The Desbrum et al.'s dis- 
cretization (jlOl) has been investigated under some very restricted conditions. It 
is shown in [lOJ [121 [10] that the discretization (|10l) converges to the LB oper- 
ator under the conditions that the valence of the vertex w is 6 and v = F{q), 
Vi — F{qi) for a smooth parametric surface F and the relations (7^+3 + Qi = 2q, 
? = 1, 2, 3 hold, where g^, z = 1, 2, • • • , 6 are one-ring neighboring vertices of Qi 
in the 2D domain. See [TO] [HJ [20] for more details. 



3 A new convergent discrete algorithm for LB 
operators 

In this section we will describe a simple and effective method to define the 
discrete LB operator on functions on a triangular mesh. The primary ideas 
were developed in Chen, Chi and Wu[ll [S] where we try to estimate discrete 
partial derivatives of functions on 2D scattered data points. Indeed, the method 
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that we use to develop our algorithm is divided into two main steps: First, we 
lift the 1-neighborhood points to the tangent space and obtain a local tangential 
polygon. Second, we use some geometric ideas to lift functions to the tangent 
space. We call this a local tangential lifting (LTL) method. Then we present 
a new algorithm to compute their Laplacians in the 2D tangent space. This 
means that the LTL process allows us to reduce 2D curved surface problems to 
2D Euclidean problems. As one will see later, our approach of discretization is 
quite different from the discretizations discussed in section 2. 

Consider a triangular surface mesh S = {V,F), where V = {vi\l <i< ny} 
is the list of vertices and F = {//cjl < k < np} is the list of triangles. 



3.1 The local tangential lifting (LTL) method 

To describe the local tangential lifting (LTL) method, we introduce the approx- 
imating tangent plane TSa{v) and the local tangential polygon Pa{v) at the 
vertex v of as follows: 

1. The normal vector Na(v) at the vertex w in S* is given by 

«,(„) = „^"^'-' (12) 

II T,T€T{v) Wt^VtII 

where T{v) is the set of triangles that contain the vertex v, Nt is the unit 
normal to a triangle face T and the centroid weight is given in [51 [3] by 

1 

i^T = ^ 1 (13) 

where Gt is the centroid of the triangle face Tdetermined by 

Gr = (14) 

Note that the letter A in the notation Na{v) stands for the word "Ap- 
proximation" . 

2. The approximating tangent plane TSa{v) of S* at w is now determined by 
TSa{v) = {wG R^\w±NAiv)}. 

3. The local tangential polygon Pa{v) of v in TSa{v) is formed by the vertices 
Vi which is the lifting vertex of Vi adjacent to v in V: 

v, = {v,-v)-<Vi-v,NA{v)>NA{v) (15) 

as in figured] 

4. We can choose an orthonormal basis ei,e2 for the tangent plane TSa{v) 
of 5' at w and obtain an orthonormal coordinates {x, y) for vectors w G 
TSAiv) hy w — xei + ye2- We set Vi = XiCi -\- 2/^62 with respect to the 
orthonormal basis ei, 62- 
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^ ^V,(v)' 










Local tangential 
polygon (v) 



Figure 2: The local tangential polygon Pa(v) 



Next we explain how to lift locally a function defined on V to the local 
tangential polygon Pa{v). Consider a function h on V. We will lift locally the 
function h to a function of two variables , denoted by h, on the vertices Vi in 
PAiv) by simply setting 

h{xi,yi) ^ h{v.i) (16) 

and /i(0) = h{v) where is the origin of rS'^(t;) . Then one can extend the 
function /i to be a piecewise linear function on the whole polygon Pa{v) in a 
natural and obvious way. 



3.2 A new discrete 2D Laplacian algorithm and configu- 
ration equation 

In this section we present a new discrete 2D algorithm for Laplacians acting on 
functions on the 2D domains in the x — y plane. Given a function / on a 
domain 17 in the x — y plane with the origin (0, 0) G 12, Taylor's expansion for 
two variables x and y gives 

2 2 

/(x, y) = /(O, 0)+a:/,(0, 0)+y/,(0, 0)+y /,,(0, 0)+xy/,,(0, 0)+|-/,,(0, 0)+O(r3) 

(17) 

when r = a;^ + j/^ is small. 

Consider a family of neighboring points (a;^, yj) G 17. i = 1, 2, • • • , n, of the 

n 

origin (0,0). Take some constants ai, i — 1,2,- •• ,n, with of — Then 
one has 

n 

Ea,(/(a:„ 2/0-/(0,0)) 
j=i 

= (i:a.x,0/.(0,0) + (Ea.2/0/y(0,0) + i(Ea,a;f)/,,(0,0) (18) 

2=1 2—1 i—1 

n n 

+ (E a,x,y,)/.y(0,0) + i(E a,y2)j^^(o^ q) + ©(r^) 

4=1 1=1 

We choose the constants a^, i = 1, 2, • • • , n, so that they satisfy the following 
equations: 



6 



aiXi = 0, 



^a^y^ = 0, 



(I) 



(11) 



and 



aiXiVi = 0, 



(III) 



4=1 



or equivalently 



X? - yf) = 0, 



(IV) 



One can rewrite these equations in a matrix form and obtain the foUowing 
the configuration equation: 



2/1, 



X2, 

y2, 



Vn 



Xiyi, X2y2, ■ ■ • , Xnyn 

\Xi ~ Uii X2 — y2: ■ ■ ■ , Xj^ — y^J 



a2 



\oin/ 








(19) 



The solutions ai of this equation allow us to obtain a formula for the Laplacian 
A/(0,0): 

A/(0,0) =/..(0,0) + /yy(0,0) 



2 E «.(/(^.:y.)-/(0.0)) 



For the reason of symmetry, Equation ([^(1)) also gives 



0(r) 



(20) 



4Ea,(/(x„y,)-/(0,0)) 

A/(0, 0) = + 0(r) 

E ai{xl + yj) 



(21) 



since we have ^ ctixf = ^ cti2/^ • It is worth to point out that Equation (j2ip 

i=l 1=1 

is an generalization of the well-known 5-point Laplacian formula. In the 5-point 
Laplacian case, we have the origin (0, 0) along with 4 neighboring points (s, 0), 
(0, s) (— s, 0) and (0, —s) for sufficiently small positive number s. From Equation 
([TO)) , one can find a solution a; = ^ for i = 1, 2, 3, 4, in this case. 



7 



3.3 A new discrete approximation for LB operators over 
surfaces 

Now we can come back to handle the local lifting function h and propose a new 
discrete approximation for the LB operator over a triangular surface mesh S . 
From Equation (l^lT) . we can define a discrete LB operator for the function 
h at the vertex v by 

n 

AaHv) = (22) 

i=l 

where the constants a^, i ~ 1,2,-- - ,n satisfy the configuration equation (|19p. 
Note again we have Vi = [vi — v)— < Vi — v, Na{v) > Na{v) and Vi — 0:^61+2/^62. 
Indeed, this definition of the discrete LB operator is independent of the 
choice of the orthonormal basis 61, 62- It depends on the choice of the constants 
ai. To obtain a unique solution for i — 1,2,-- - ,n, with Q^l — 1j oi^^ 

can simply choose 5 closest neighboring points Vi of the origin in Pa{v). In case 
that the local polygon has less than 5 vertices, one can also lift the 2 or 3-ring 
of neighboring vertices of v in V. In this way, we can call A^/i(w) a 6-point 
Laplacian formula. 

3.4 A convergence theorem for the discrete LB operators 

In this section we will prove that the discrete LB operator A Ah for a smooth 
function /i on a regular surface E is a convergent 0{r) approximation of the true 
LB operator A-^h. To show this, let E be a smooth regular surface in the 3D 
Euclidean space and p G E. Consider the exponential map expp : rE(p) E 
from the tangent plane TYi{p) of E at the point p into the surface E. See do 
Carmo [8l[9] for discussions about the properties of the exponential map expp. 
One of the well-known properties of the exponential map exp^ is that it is a 
local diffeomorphism around the origin G TE(p). In other words, if W is 
a sufficiently small open domain around the origin 0, expp : W ^ D \s a 
diffeomorphism where D = exp(M^) is an open domain around p. In particular, 
the inverse of exp^ exists on D. 

Given a smooth function h on the regular surface E, we can lift h via the 
exponential map exp^ locally to obtain a smooth function h defined on g 
TE(p) by setting 

h{w) = /i(expp(w)) (23) 

for w € W. Fix an orthonormal basis 61,62 for the tangent space TE(p). 
This gives us a coordinate system on TE(p). Namely, for w € W we have 
w = xei + ye2 for two constants x and y. Without ambiguity, we can identify 
the vector w & W with the vector {x,y) with respect to the orthonormal basis 
61,62. In this way, the function h can also give us a smooth function h of two 
variables x and y by defining 

h[x,y)^h{w) (24) 
for w = xei + ye2- Using these notations, we will prove 
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Lemma 1. One has 

A^hip) = Ah{0) = Ah{0, 0) (25) 

Proof: 

It is well-known that the LB operator A^h{p) acting on a smooth function 
at a point p can be computed from the second derivatives of h along any two 
perpendicular geodesies with unit speed. See do Carmo [5] for details. Indeed, 
we consider the following two perpendicular geodesies with unit speed in E by 
using the orthonormal vectors 61,62' 

6i{t) ^ expJtei), 1 = 1,2 (26) 



with Ci{0) = p and ^(0) = e,;. One has 



=A^(0) (27) 
=A/i(0,0) 

□ 



Next we consider a triangular surface mesh S = {V, F} for the regular surface 
E, where V = {vi\l < i < ny} is the list of vertices and F = {Tk 1 < k < np} 
is the list of triangles and the mesh size is less than r. Fix a vertex v in V. For 
each face T Cz F containing v, we have 

N^{v) ^ Nt + 0{r^) (28) 

where N^{v) is the unit normal vector of the true tangent plane TT,{v) of E at w 
and Nt is the unit normal vector of the face T. Since the approximating normal 
vector Na{v), defined in section 3.1 is a weighted sum of these neighboring face 
normals Nt, we have 

Lemma 2. One has 

N^{v) = Nav + Oir^) (29) 

Due to this lemma, the orthonormal basis ei, 62 for the tangent plane TE(w) 
will give us an orthonormal basis 61,62 for the approximating tangent space 
T^Aiv) = {w e R^\wJ-NAiv)} by the Gram-Schmidt process in linear algebra: 

^ ei- < ei,NAiv) > Na{v) 
\\ei- <ei,NA{v)> NA{v)f 

and 

_ 62" < e2,NA{v) > Na{v)~ < 62,61 > ei 
~ \\e2- < e2,NA{v) > Na{v)- < 62, ei > 6i|| ' 
Logically speaking, one can first choose an orthonormal basis 61,62 for the 
approximating tangent space TS'^(w) and then apply the Gram-Schmidt process 
to obtain an orthonormal basis ei, 62 for the tangent plane TE(w). In either way, 
we always have by Lemma [2] the following relations. 
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Lemma 3. One has 

e, =e, + 0(r2), i = l,2 (30) 

Consider a neighboring vertex Vi of v in V . For r small enough, we can use 
the inverse of the exponential map exp^ to lift the vertex Vi up to the tangent 
plane rS(i;) and obtain 

ill = exp;;^(wi) e ty,[v) 

and 

Vi = x^ei + yie2 

for some constants. As discussed in section 3.1, we can also lift the vertex Vi up 
to the approximating tangent plane TY,a{v) and get 

Vi = [vi - v)- <Vi-v, NAiv) > NAiv) 

and 

Vi = x^ei + yie2 
for some constants Xi,yi. Then Lemmas [2] and [3] yield 
Lemma 4. One has 



(31) 



j ii =x^ + Oir^) 

Using these relations, one can solve the configuration equation (jl9p for 
(ii, yi) and {xi, yi) respectively and obtain their corresponding solutions oii and 
Ui with the relation 

5, =a, + 0(r2) (32) 

Note that the lifting function h is a smooth function of two variables x and 
y. Equation (j2ip in section 3.2 now gives an approximation of the Laplacian 
A/i(0,0): 

n 

AY.a,{h{x^,y,)-h{Q,Q)) 

Ah{Q,0)^^^ + 0{r) (33) 

aiixi"^ + yi"^) 

i=l 

The relations ^ anA^ imply 

n 

AY. a,{h{v,) - h{v)) 

Ah{0,Q)^^ + 0{r). (34) 

E oiiix^ + yf) 
1=1 

This along with Lemma [T] proves the following convergence theorem. 

Main Theorem. Given a smooth function h on a regular surface S and a 
triangular surface mesh S — (V, F) with mesh size r, one has 

A^h{v) = AAh{v) + 0{r) (35) 

where the discrete LB operator AAh{v) is defined by Equation 



10 



Remark 1. The discussions in this section also indicate that as long as we 
have 0{r) -convergent algorithms to estimate gradients, Laplacians and other 
intrinsic derivatives of 2D smooth functions, the LTL method and methods in 
section 3.4 will allow us to develop corresponding discrete convergent algorithms 
over 3D surfaces. It is possible to obtain a 0(r^) algorithm by extending Taylor's 
expansion {IT^to the third order and improving the configuration equation hl9\l . 



4 Numerical simulations 

In this section, we shall compare two convergent Laplace-Beltrami methods: 
Xu's method (see [2^) and our proposed method. We take four functions, 

Fi{x,y) = v/4- (0:^0.5)2 -(2/ -0.5)2. 



F2{x,y) = tanh(9x — Oy). 

p , X _ 1.25+cos(5.4jy) 
^3[-^,y) — 6+6(3x-l)2 ■ 



F,{x,y) = exp {-%{{x - 0.5)2 + _ 5)2)^ 
over xy-plane as three dimensional surfaces. 




Figure 3: The triangulation of the domain. 



The exact and approximated mean curvatures are computed as some selected 
interior domain points {xi,yj) with xi^yj ^ {0,1}. The domain (a) is a three 
directional triangular partition, the domain (b) is a four directional triangular 
partition and the domain (c) is a unstructured triangular partition. To observe 
the convergence property, these domains are recursively subdivided by the bisec- 
tion linear subdivisions. Hence, h = ~, i = 1, 2, • • • , where Hq — -\/0.2, 0.23 
are the maximal value of edge lengths of the triangulations (a), (b) and (c), 
respectively. 

The maximal errors of these simulations are shown in the Table [T] Table 
[2] shows the time costs for the computations in the domain (c) with h = ^ . 
Obviously, our proposed method is more accurate and faster than Xu's method. 
Furthermore, the convergent rate of our method is also better than Xu's method. 
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Table 1: The maximal errors of Laplacian 



Xu's method Our method 





Fi 


2.ME - 03 * /i^ 


2.36£; - 


05*/i4 


Domain (a) 


F2 


1.81£; + 01*/i2 


3.15E;- 


02 




F3 




1.58E- 


02*/i3 




Fi 


1.50i;-+l*/i2 


1.68£; + 


00 




Fi 


4.24^; - 03 * /i^ 


2.64^; - 


05 *h'^ 


Domain (b) 


F2 




5.13i;- 


01*/i2 




F3 


1.47E + 01*/i2 


4.29£; - 


01 * h"^ 




Fi 


3.51£; + 01 * 


7ME - 


01 




Fi 


1.49i;-01 


1.53E- 


03*/i^ 


Domain (c) 


F2 


l.ME + 00*hi 


1.59E- 


01 * 




Fi 


9.22E - 01 


3.00E - 


01 * h 




Fl 


r).41£' - 01 * /) 


1ME + 


01 * //- 



Table 2: Time costs for the computations of domain (c) 



(seconds) 


Xu's method 


Our method 


Fl 


0.024 


0.010 


F2 


0.026 


0.014 


Fs 


0.026 


0.015 


Fi 


0.025 


0.012 
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